The role of final state correlation in double ionization of helium: a master equation approach 
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The process of nonsequential two-photon double ionization of helium is studied by two complementary nu- 
merical approaches. First, the time-dependent Schrodinger equation is solved and the final wave function is 
analyzed in terms of projection onto eigenstates of the uncorrelated Hamiltonian, i.e., with no electron-electron 
interaction included in the final states. Then, the double ionization probability is found by means of a recently 
developed approach in which the concept of absorbing boundaries has been generalized to apply to systems 
consisting of more than one particle. This generalization is achieved through the Lindblad equation. A model 
of reduced dimensionality, which describes the process at a qualitative level, has been used. The agreement be- 
tween the methods provides a strong indication that procedures using projections onto uncorrelated continuum 
states are adequate when extracting total cross sections for the direct double ionization process. 

PACS numbers: 32.80.Fb, 32.80.Rm, 31.15.-p, 02.70.-c 
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I. INTRODUCTION 

The problem of multiphoton single- and multiple ionization 
of atoms and molecules has been subject of intense research 
in recent years. The development of new high-frequency light 
sources, such as free-electron lasers (FEL) and high-order har- 
monic generation sources, capable of generating intense co- 
herent radiation, is one reason for this. Pulses of attosec- 
ond durations are now routinely produced by high-order har- 
monic generation with intensities as high as lO'"* W/cm^ |[l|], 
enabling experimental studies of nonlinear phenomena in the 
extreme-ultraviolet regime in atoms |0, Ht] and molecules 101 • 
Owing to recent advances in FEL technology, femtosecond 
pulses of unprecedentedly high intensity are now available, 
covering wavelengths ranging from vacuum ultraviolet to soft 
x-rays JslSl- These pulses have opened the door for experi- 
mental studies of multiphoton multiple ionization of complex 
atoms and atom clusters jT^ . A parallel development of 
ab initio numerical methods, capable of addressing the prob- 
lem of intense-field multiphoton ionization in one- and two- 
electron systems has taken place ifTll - fTsIl . 

Despite these developments, the problem of nonsequential 
(direct) two-photon double ionization of helium still remains 
an unsolved problem, even though it has been widely stud- 
ied during the last decade both theoretically lfl4l - l28ll and ex- 
perimentally [3, 29-31]. The main obstacle in experiments 
has been to produce sufficiently high ionization yields. On 
the theoretical side, accurate predictions for the total (gener- 
alized) cross sections of the process remain elusive, as values 
obtained with different methods vary by more than an order 
of magnitude. What characterizes this particular two-photon 
process is that it depends upon exchange of energy between 
the outgoing electrons, i.e., it is a so-called nonsequential or 
direct process, as opposed to a sequential ionization process 
where both electrons may be considered as independent par- 
ticles. The sequential process is energetically inaccessible for 
photon energies below 54.4 eV, but becomes the dominant 
two-photon ionization process at higher energies. Concerning 
the nonsequential process, the great discrepancies that remain 



between different theoretical calculations are usually ascribed 
to the different ways electron correlations are handled. It has 
furthermore been claimed, by several authors, that the two- 
photon double ionization process in helium is a problem in 
which electron correlations in the final state play an extremely 
important role. This stands somewhat in contrast to the related 
(direct) one-photon double ionization process, where a com- 
plete agreement between different theoretical approaches and 
experiments is now established |19, 32-34;]. 

As it turns out, it is especially the separation of the two- 
photon single ionization, where the remaining electron is left 
in an excited ionic state, from the two-photon double ion- 
ization that is the bottleneck, the main problem being that 
electrons are emitted with similar energies in the two pro- 
cesses. Moreover, the fact that the single ionization process 
is much more probable makes the problem even more chal- 
lenging, and an exact knowledge of the role of correlation in 
the final state particularly important. The effect of electron 
correlations in the final state on the total integrated cross sec- 
tion was studied in detail by Foumouo et at. I19L i26l] . Us- 
ing the J-matrix method they were able to incorporate corre- 
lation effects in the final (single continuum) scattering states, 
with the result that their calculated generalized cross section 
increased significantly compared to the result they obtained 
when projecting the final wave packet directly onto uncorre- 
lated products of Coulomb waves (with nuclear charge Z = 2) 
(compare full green curve with diamonds with full black curve 
with crosses in Fig. [Til. More importantly, their uncorrelated 
result is in complete accordance with calculations performed 
by others SIllMllSlllllllli, which all have in com- 
mon that the role of correlation was completely neglected in 
the final state. For comparison, these results, together with the 
result of Feist et al. ll24fl and Nepstad et al. ll28ll . and available 
experimental data, is shown in Fig. [1] 

Generating correlated multichannel states, and employing 
lowest-order (non-vanishing) perturbation theory, Nikolopou- 
los and Lambropoulos [23] obtained a value for the cross sec- 
tion that differed even more substantially from the uncorre- 
lated results (black line with triangles in Fig. [T]!. From this 
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FIG. 1 : (Color online) Two-photon double ionization cross sections. 
Black circle: experimental result of Hasegawa et al. l29h . red cross: 
experimental result of Sorokin et al. 13011 . blue line with circles: the 
results obtained by Nepstad et al. I2al . red line with squares: Feist 
et al. l24ll . green line with diamonds: Foumouo et al. jT^Kwith cor- 
relation, WC), black line with crosses: Foumouo et al. 11911 (no cor- 
relation, NC), orange dashed line with '+' (scaled): ID result with 
the Hamiltonian and black line with triangles: Nikolopoulos et 
al. ll23ll . The vertical lines define the two-photon direct double ion- 
ization region. 



observation, it appears reasonable to conclude that the inclu- 
sion of electron correlation in the final states is required. How- 
ever, the R-matrix calculations by Feng and van der Hart ifTsll . 
calculations employing exterior complex scaling by Horner et 
al. ll22ll , and the time-dependent close-coupling calculations 
by Ivanov and Kheifets lulll , which all consider final state cor- 
relations, produce results that are similar to the uncorrected 
results. Thus, at present, the matter remains unresolved. 

In this article, we revisit the problem of nonsequential two- 
photon double ionization of helium. We employ an approach 
based on the concept of absorbing boundaries and the Lind- 
blad formalism ||35I1 to calculate the generalized cross sec- 
tion for the process. Absorbing boundary conditions are often 
quite useful in the context of wave packet propagation ll36l - 
13811 . However, applying them in the study of systems consist- 
ing of several particles is not straightforward. This is due to 
the fact that when one particle is absorbed, the Schrodinger 
equation does not retain any information about the remain- 
ing particles. Here we will demonstrate how the notion of 
absorbing boundary conditions may be generalized to apply 
to many-particle systems in such a way that single and double 
ionization may be distinguished Issll . Specifically, the method 
allows for describing the dynamics of the remaining electron 
as the first one is absorbed due to photoionization. A distinct 
advantage is that an exact knowledge of the final scattering 
states is superfluous. Furthermore, a direct comparison with 
the result obtained by projecting the final wave function onto 



eigenstates of the uncorrected Hamiltonian reveals that this 
procedure indeed is adequate, as a complete agreement be- 
tween the two different methods is demonstrated. 

Describing a helium atom exposed to a laser pulse numer- 
ically is a challenging task both in terms of algorithm and 
computational power. Several theoretical studies of this sys- 
tem involves models of reduced dimensionality ll39l - l43ll . Only 
during recent years has a full three-dimensional (3D) descrip- 
tion become numerically feasible. The formalism based on 
the Lindblad equation involves resolving the dynamics of a 
one-particle density matrix in addition to a two-particle wave 
function - both involving six degrees of freedom when de- 
scribing the full 3D problem. Since these two sub-systems 
are coupled, more involved numerical schemes are necessary 
when solving these equations than in the case of the time- 
dependent Schrodinger equation (TDSE). Thus, we resort to 
a one-dimensional (ID) model atom in the present work. We 
argue, however, that the coincidence between the double ion- 
ization probabilities obtained by the two different methods do 
shed light on the question of whether projection onto uncorre- 
cted final states is an adequate approach also in the 3D case. 

Having validated the Coulomb-wave projection method, we 
next consider the convergence property of the cross section 
in the vicinity of the threshold at 54.4 eV. It has been re- 
ported that the cross section exhibits an apparent divergence 
in this Hmit S HHH HHil, a behavior that is usually 
interpreted as being due to an unwanted inclusion of the se- 
quential process in the calculations caused by the bandwidth 
of the pulse ll45l l46ll . However, it has recently been shown 
that this rise cannot solely be attributed to such an effect 12^ . 
That 3D study was limited to photon energies below 52 eV. 
Going beyond the 52 eV limit would require the application 
of extremely long pulses with durations of several tens of fem- 
toseconds, which ultimately becomes an unmanageable prob- 
lem in 3D. A ID model remains computationally tractable, 
however, even for these long pulses, thus allowing to study 
the behavior of the cross section very close to the threshold. 

In the following section the theory is outlined. Specifically, 
Sec. Ill Al describes the model and discusses to what extent 
it is able to make predictions also valid for the full 3D sys- 
tem. The two different methods used in order to obtain to- 
tal cross sections are described in Sees. IIIBI and III CI Their 
implementation, and possible extension to 3D is discussed in 
Sec. HID I The comparison between the two methods is pre- 
sented in Sec. |III] Furthermore, the behaviour of the cross 
section when the photon energy approaches the threshold at 
54.4 eV is discussed. Conclusions are drawn in Sec. lIVI 



II. THEORY 
A. The model atom 

The Hamiltonian of the system is given by 

H = /l(A-l) + /2(x2)+yint(xi,X2) (1) 
1 

h{x) = --^+yNuciW+.«£(f) (2) 
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where x,- is the position of particle / and E is some external 
time-dependent electric field of finite duration T . Here, we 
have introduced atomic units, denoted "a.u.", which are de- 
fined by choosing the electron mass, the electron charge and 
h as units for their respective quantities. Both the interaction 
with the nucleus, VnucI, and the electron-electron interaction, 
Vint, are taken as regularized Coulomb potentials, instead of 
bare Coulomb potentials, i.e.. 



Vnuci(-«) = - 



^Nucl 



Zint 



(3) 
(4) 



By choosing ZnucI = 1.1225, Zjnt = 0.6317 and 5 ~ 
0.3028 a.u., the ground state and the first and second ioniza- 
tion thresholds coincide with those of the 3D case. We em- 
phasize that it is crucial that also the second ionization thresh- 
old, i.e., the first excited state of He+, is correctly reproduced. 
For the photon energies considered here, it means that one- 
photon excitations of He+ is not possible, concordant with the 
3D case; indeed, such a possibility would introduce artifacts 
detrimental to our present purpose. 

While these are necessary conditions for the model to pro- 
duce results comparable with the original system, we cannot 
take for granted that they allow for conclusions pertinent to 
the real 3D to be drawn. Crucial geometrical aspects may 
potentially be lost when reducing the dimensionality of the 
problem. For instance, in the case of one-photon double ion- 
ization of helium, it is known that back-to-back scattering of 
the two electrons is prohibited for equal energy sharing ll26ll . 
Moreover, electron ejection is more probable in directions per- 
pendicular to the polarization axis than parallel to it. Thus, for 
such a process, our ID model cannot be expected to provide 
relevant information. For two-photon double ionization, how- 
ever, it is known that back-to-back emission along the polar- 
ization axis is the dominating mode of ionization 12611 . Such a 
process is indeed described within the ID model. 

Predictions for the angular aspects of the double ionization 
process, such as the directional distribution of the photoelec- 
trons, cannot be made within the model. However, when it 
comes to the total cross section, it has been argued that for 
the process at hand, radial correlations, as opposed to angular 
correlations, provide the dominant contribution ||47ll . Hence, 
it seems reasonable to expect that the model should be able 
to describe the non-sequential double ionization process at a 
qualitative level. This assumption is strengthened by the fact 
that the ID cross section, albeit too low in absolute value, does 
resemble the full 3D cross section as far as the photon energy 
is concerned, see Fig.[T] 

In addressing the issue of whether projection onto uncorre- 
cted final states provide the correct double ionization prob- 
ability asymptotically, it is crucial not to underestimate the 
role of the electron-electron interaction. In the asymptotic re- 
gion, both the nuclear potential and the electron-electron re- 
pulsion of the model coincide with the 3D case. Moreover, 
ID models actually tend to overestimate correlation. This 
can be understood from the fact that the "smoothing" of the 



Coulomb potential, Eq. (O, effectively reduces the kinetic en- 
ergy of the system, which in turn increases the importance 
of correlation ll48ll . In our model, this is confirmed by com- 
paring the expectation value of the kinetic energy T for the 
ground state with that of the original system; in ID we find 
that (r) = 0.685 a.u., which is about a fourth of the value 
found in 3D. The same tendency is found when comparing 
the correlation energy of the ground state; a Hartree-Fock cal- 
culation for the ID model gives a ground state energy which 
is off by about 18 % as opposed to 1.4 % in the 3D case ll49ll . 
Thus, we expect electron-electron correlation to be at least as 
important in the ID model as in 3D. 

We note that, as the Hamiltonian is independent of electron 
spin, and the initial state is a spin eigenstate with symmetric 
spatial part, the system may be described formally as a two- 
boson system without spin. This will simplify the notation in 
the following. 



B. Uncorrelated final state approach 



In the method of proj ection onto uncorrelated final states, 
following Refs. 01 [iMl US IM IH IH, we start out by 
diagonalizing the one-particle Hamiltonian /; of Eq. (|2]l, h (p„ = 
e„^„. The "uncorrelated double continuum" is defined by the 
span of two-particle symmetric states with positive energies, 

^^^7^ ^N,„,„[(j)„,(xi)(j)„{x2) + (j)„{xi)(j),„{x2)] (5) 



with £,„ , £„ > and N„, „ = < \{2 ' 



m <n 



Clearly, these states do not represent the actual double contin- 
uum states of the system. However, assuming that the electron 
correlation diminishes in significance as the doubly ionized 
wave packet travels outwards, the correct double ionization 
probability should be obtained by projecting the unbound part 
of the final state wave function at sufficiently long time after 
the interaction onto the uncorrelated continuum states. 



(6) 



m n>m 



where the projection operator removes the bound part. 



C. Absorbing boundaries and the Lindblad equation 

We seek to determine the validity of the above approach 
by comparing the double ionization yields thus obtained to 
the ones obtained using absorbing boundaries. These are in- 
troduced by augmenting the Hamiltonian with a complex ab- 
sorbing potential, — ;T(x). The function r(x) is zero within 
a certain region of x and beyond this region it is positive and 
increasing with distance. 



Y{x) = 0, 
T{x) > 0, 



\X\ <XT 
\x\ >XT. 



(7) 
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Waves entering into the region where F > are attenuated and 
eventually die out completely. 

Simply introducing the absorber in the TDSE is problem- 
atic when we wish to distinguish between single and double 
ionization. This becomes evident when considering the fact 
that the wave function is normalized to the probability of all 
particles being represented. If one electron travels outwards 
and is subsequently absorbed, not only is all information about 
this electron lost, but so is all information about the remain- 
ing electrons as well. Thus, only information about total ion- 
ization probabilities may be retained when using an absorber 
in combination with the TDSE, while distinguishing between 
single, double etc. ionization is hard to achieve. 

In a recent paper it was demonstrated how the notion of ab- 
sorbing boundaries may be generalized to an A^-particle con- 
text in such a way that the remaining (A^— 1) -particle system 
is recovered after one particle has been absorbed from the sys- 
tem Issll . Likewise, the (A^ — 2) -particle system is recovered 
when yet another particle is absorbed etc. It was argued that, 
since this is a Markovian process, the Lindblad equation is the 
proper framework within which to achieve this ifsollsill . In 
general, the dynamics of a system where particles are lost due 
to a complex absorbing potential is described by Issll 



ih—pn = [i/,p„]-/{f,p„} 
at 



2//r(^)v/(^)p„+iV/+(^)fl'^, 



(8) 



where p„ is the density operator corresponding to the n- 
particle sub-system, « ~ 0,1,.. .,N. The generalized coordi- 
nates £, correspond to all degrees of freedom. The operators 
and V'^^(^), which annihilate and create, respectively, 
a particle with position and spin given by ^ , obey the usual 
anti-commutation relations for fermions: 

{V/(^),V/(r)}-0, {v/(^),v/^(r)} = 5(^-r)- (9) 

The operators H and F are the Hamiltonian and the complex 
absorbing potential, respectively, expressed in terms of these 
field operators such that their explicit forms do not depend on 
the particular number of particles. Specifically, for the two- 
particle system at hand the evolution may be written as 



in-\'¥2) = (//-/r(xi)-/r(x2))|i'2) (lO) 

at 



at 



2i rix)xix)\^2){'V2\xHx)dx 



^^Po = 2 / r{x)x{x)pixHx)dx, 



(11) 



(12) 



where |^2) is the two-particle state, pi is a one-particle 
density operator and po = Po{t)\^){—\ is a density opera- 
tor corresponding to vacuum, i.e., no particles present. In 
Eqs. ( I10I11I12T | the spin degree of freedom is integrated out. 
Hence the field operators j'^' (x), which annihilate (create) a 
particle in position x, correspond formally to spinless bosons. 



We emphasize that this formalism does not rely on the reduced 
dimensionality of the model. The equations remain the same 
regardless of the dimensionality of the coordinates x (or ^). 

We see that I't follows the TDSE. Equation ( fTTT i. which 
governs the evolution of the one-particle sub-system, has sev- 
eral terms. The first term on the right hand side, which cor- 
responds to the von Neumann equation, simply governs the 
evolution of pi dictated by the one-particle Hamiltonian h. 
The second term corresponds to removal of the particle due 
to the absorber, whereas the last term is a source term de- 
pending on the overlap between the absorber and The 
decrease in norm of 'Pt is identical to the increase in the trace 
of pi induced by this source term, and the remaining electron 
is properly reconstructed. Similarly, the vacuum-probability 
PQ increases by the same amount as Tr(pi) decreases due to 
absorption. These mechanisms are illustrated schematically 
in Fig. |2] The sum of all norms and traces remains equal to 
unity at all times ifsill . i.e., 

P2{t) + piit)+po{t) = l, Vf (13) 
with p2{t) = \\^'2{t)f and pi(f) = Tr(pi(f)). 

As also illustrated in Fig. |2l the density operator pi de- 
scribes in general a mixed state. 



(14) 



i.e., one cannot in general expect a single-particle wave func- 
tion since correlation information is lost whenever particles 
are removed. 

In Fig.[3]we have shown p2{t), pi (f ) and po{t) as functions 
of time both during and after the interaction with the laser 
pulse. It is seen that the two-particle probability decreases 
during and after the pulse, whereas the one- and zero-particle 
probabilities increase - all towards some asymptotic value. It 
is also seen that po{t) converges much more slowly than P2{t) 
and pi{t). We will return to this issue later. In the upper panel 
the probability for the system to remain in the ground state 
is also displayed. It shows that practically all population of 
bound states is confined to the ground state; very little excita- 
tion is seen. 

The key point here is that as the first electron is absorbed, 
the dynamics of the remaining one may still be described. 
Moreover, if also the second electron is absorbed, we will 
eventually find the total system in the vacuum state. Now, sup- 
pose we may choose the absorption-free region large enough 
for the probabilities po, pi and p2 to converge towards val- 
ues independent of xj-, then these probabilities are subject to 
straightforward interpretations: p2{t °°) is the probability 
that the system remains bound after interaction, pi{t ^ °o) is 
the single ionization probability, and p(j{t — > oo) is the double 
ionization probability. This interpretation relies on the idea 
that there is a finite distance from the nucleus beyond which 
we may conclude that the electron is indeed in the continuum. 
The validity of this assumption, in turn, rests on whether the 
probabilities pi{t — > °o) and po(t °o) converge as the ab- 
sorber is moved outwards. 

The ideal absorber should be both transmission and reflec- 
tion free, and hence not influence the dynamics in the region 
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FIG. 2: (Color online) A schematic view of the absorbing bound- 
ary multi-particle method. The system is initially described by a 
two-particle wave function (top). As a particle is ionized and subse- 
quently hits the absorber near the boundary (F), the remaining elec- 
tron is described by a one-particle density matrix pi (middle level). 
Some of the eigenstates of pi are shown. Further ionization cause 
also pi to overlap with the absorber, which in turn causes the vac- 
uum state to be populated. See text for further details. 



where r(x) = at all. In a many-body context, a certain indi- 
rect influence is difficult to avoid, however, since removal of 
one particle also removes the interaction between this particle 
and the remaining ones. This may be crucial if the field moves 
the electrons such that they are close while one of them is be- 
ing absorbed. However, given a large enough xj and long 
enough propagation time, the electron-electron repulsion it- 
self should cause the importance of this interaction to dimin- 
ish. Hence we conclude, once again, that the approach should 
be valid if one is able to obtain absorber-independent results. 

As in the case of projection onto uncorrected final states, 
the state of the system must be propagated a certain time be- 
yond the duration T of the electric pulse in order for the dou- 
ble ionization probability pQ to reach convergence. Having 
obtained the probability of direct double ionization, Pqi, for a 
situation in which perturbation theory applies, the total cross 
section for direct double ionization is found as ll24ll52ll 



— — , ^eff = 



lit) 
k 



dt. 



(15) 



where /(f) is the intensity of the laser pulse, and Iq is the max- 
imum intensity. For a pulse with a sine square envelope. 
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FTG. 3: (Color online) The probability of having two (upper panel), 
one (middle panel) or zero particles (lower panel) on the grid. The 
thin, black curve in the upper panel shows the probability for the 
system to remain in the ground state. The pulse duration is here 
r = 118 a.u. = 2.86 fs, the peak intensity is 2.2 x lO'^W/cm^, 
and the central frequency of the pulse corresponds to a photon 
energy of 43.5 eV. We see that both the one-electron probability, 
p\{t) =Tr(pi(/)), and the zero-particle probability po{t) increase 
as the norm of the two-particle probability, P2(0 = ||1'2(0II^' de- 
creases. The dashed curve in the lower panel represents the double 
ionization probability obtained by analyzing the absorbed part of the 
wave function "on the fly". See text for details. 



E{t) ^Eosin^ 



Ttt 



sm{m + (p), 



(16) 
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we have T^g = 357/128. 



D. Numerical considerations 

As mentioned in the introduction, one of the reasons 
for choosing a ID model atom rather than trying to solve 
the full 3D problem is related to the implementation of 
Eqs. ( I10I11I12] |. While there are no formal obstacles prohibit- 
ing the application of the method to a 3D situation, and in 
spite of the fact that a grid of relatively small extension may 
be used, such calculations would generally involve very large 
basis sets. Thus, in a practical situation a parallelized version 
of the scheme is required. Since the two-particle dynamics 
is simply dictated by the Schrodinger equation, this part may 
be handled by any of a number of proposed parallel numeri- 
cal schemes available in the literature |24, 28, 53*], given the 
availability of sufficient computational resources. A similar 
scheme could then be adapted for the one-particle density ma- 
trix, Eq. ( fTTT i. The remaining bottleneck is the calculation of 
the source term, which must be updated at each time step. Ev- 
ery part of the two-particle wavefunction at the previous time 
step which has a finite overlap with the absorber will con- 
tribute to the source term. This operation is potentially quite 
demanding on communication between the processing ele- 
ments (PE), and some scheme to partially assemble the source 
term with the local data on each node is probably required, 
before these sub-elements may be transmitted to the correct 
PE for final assembly. Alternatively, the complexity of solv- 
ing Eq. ( fTTT ) could be reduced by a lower rank approximation, 
e.g., along the lines proposed by Leth and coworkers [Sj]. In 
any case, a non-parallel ID model implementation is a natural 
starting point for further developments of the method, produc- 
ing benchmarks and providing a convenient testing ground for 
determining optimal implementations of the various parts of 
the algorithm. 

Returning to our ID model, both the TDSE and 
Eqs. dlOll lll2T i are solved by means of split operator tech- 
niques ifssll on uniform numerical grids, xj ~ — L/2 + {j — 
I) Ax, where ; = 1,2 is the particle index, and Ax ^ L/ {,yV — 
1), where ,yV is the number of grid points. The code solving 
the TDSE, implemented in the Pyprop framework ll28ll56ll , is 
parallelized and thus able to handle relatively large grids. 

The numerical scheme applied in order to solve 
Eqs. (I10I11I12T | was presented in jssll . However, here 
we have made some adjustments to the method in order to 
accommodate for two challenges which are both related to 
the long-range nature of the Coulomb interaction with the 
nucleus. Firstly, the system features Rydberg states, which 
may have a finite overlap with the absorber. If these states are 
populated, it could induce artificial ionization on large time 
scales. Hence, at f = T all components of the bound states of 
the helium system was removed. These states, including the 
initial ground state, were found by propagation in imaginary 
time. Secondly, it is well known that photoelectrons of low 
kinetic energy may emerge from the process at hand ll26ll . 
This may be somewhat problematic since, as we have seen 
in Fig. |3] it causes po(0 to converge rather slowly in time (it 



takes the second electron a long time to reach the absorber). 
Moreover, the low energy of the photoelectron makes it 
vulnerable to artificial reflections induced by the absorbing 
potential. After the interaction with the laser pulse, t > T, 
this problem may be circumvented by analyzing the source 
term in Eq. (fTTT i "on the fly" instead of propagating pi further. 
Suppose that the "amount" of double ionization is known at a 
time t > T, then at the next time step, t + At, the source term 
in Eq. ( fTTT i adds a contribution (Api)j to the one-particle 
density operator. Since it is not very costly to diagonalize 
the one-particle Hamiltonian h, as compared to diagonalizing 
the full two-particle Hamiltonian H, this contribution (Api)j 
may be separated into a part corresponding to bound He+ and 
a (double) continuum part. The double ionization probability 
is then obtained simply by repeatedly adding (integrating) the 
latter contribution to the double ionization probability. Con- 
sequently, it is not necessary to solve Eqs. ( I11I12I ) for t > T 
- as long as the absorbed part of ^2 is properly analyzed. In 
the lower panel of Fig. [3] it is clearly demonstrated that the 
converged value for the double ionization probability is more 
easily obtained this way. 

Note that during the interaction with the laser pulse, t <T, 
it is necessary to resolve the one- and zero-particle dynam- 
ics as well, since the remaining electron is still subjected to 
the electric field. We would like to point out that the formal- 
ism based on the Lindblad equation allows for analysis of the 
ionization channel yields using a grid that is smaller than the 
actual extension of the complete ionized wave packet. This 
holds also during interaction with the laser pulse. 

For the Lindblad-formalism, practically convergent results 
for the double ionization probabilities were obtained using a 
box of size L = 130 a.u., a grid spacing of Ax = 0.254 a.u. 
and a numerical time step Af = 2.50 x 10^^ a.u. For the 
complex absorbing potential we have chosen a Manolopoulos 
form, c.f. IstIi . which has several advantages such as being 
completely transmission free and only containing one param- 
eter. We have fixed this parameter by the choosing .x: 7- = 0.2 L, 
c.f. Eq. (|7]i. Finally we note that, due to trace conservation, 
Eq. ( fT3] ). it is not strictly necessary to solve Eq. (fT2l) along 
with Eqs. ( IIOII 11 1. However, since our numerical scheme is 
not manifestly unitary, we have still done so in order to check 
that Eq. ( fT3] l indeed is fulfilled to a satisfactory degree. 

In the TDSE calculations, where the final wave packet 
is projected directly onto uncorrected products of Coulomb 
waves (with nuclear charge Z = 2), Eq. much larger boxes 
are required in order to contain the entire ionized wavepacket 
during the time evolution. We used L < 2400 a.u., which was 
sufficient for pulse durations up to 20 fs, and photon energies 
up to 53 eV. The grid resolution was fixed at Ax ~ 0.195 a.u., 
while the time step was At = 5.00 x 10^^ a.u. 



in. RESULTS AND DISCUSSION 

Figure |4] shows the total cross section for the process of 
nonsequential double ionization by two photons. Results have 
been obtained using both methods discussed above. In both 
cases, the electric field was given by Eq. (fTST l. For lower 
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FIG. 4: (Color online) The (generalized) cross section for nonse- 
quential two-photon double ionization of our ID model helium atom. 
The full curve is obtained by solving the TDSE on a large grid with 
subsequent projection of the final wave packet directly onto uncorre- 
lated products of Coulomb waves (with nuclear charge Z = 2). The 
diamonds are calculated by the method involving absorbing bound- 
aries and the Lindblad equation. 

photon energies the duration of the field corresponded to 30 
optical cycles, whereas increasingly longer pulses where nec- 
essary in order to obtain converged cross sections as the pho- 
ton energy increases towards the threshold at hco ~ 54.4 eV. 
We see that the cross sections obtained in the two different 
ways are in very good agreement. This provides support for 
the proposition that the uncorrected final states, Eq. are 
indeed able to provide accurate information about photoelec- 
trons emerging from helium asymptotically. 

In ll24ll the strong dependence of the cross section as de- 
fined in Eq. dTSI l on the pulse length T as Tico — > 54.4 eV was 
discussed. It was demonstrated that a, which seems to have 
a sharp rise, depends rather strongly on the pulse duration 
T in this limit. A possible explanation for this was consid- 
ered: for finite T, there is also a finite spectral width of the 
pulse, enabling sequential double ionization also for central 
frequencies slightly below threshold. As the probability of se- 



quential double ionization scales with T^, as opposed to 
in the nonsequential case, this gives rise to a pronounced in- 
crease in <7. The 4 fs pulse used in ll24ll . although not long 
enough to resolve the cross section behavior very close to the 
threshold, was sufficient to observe the start of an increasing 
trend, and thus rejects the sequential overlap hypothesis with 
some confidence. However, in order to resolve the true be- 
havior of the cross section in the immediate vicinity of the 
threshold, very narrow bandwidths, i.e., very long pulses, are 
required. By performing calculations with increasing pulse 
duration (up to 20 fs), we have confirmed that the value of 
a does converge towards a well-defined value as T increases 
- also near the threshold. Moreover, as the spectral width of 
the pulse becomes so narrow that the overlap with the region 
in which sequential two-photon double ionization may take 
place vanishes, the pronounced rise in the cross section is seen 
to prevail. Hence, the fact that a increases sharply as hco ap- 
proaches 54.4 eV cannot be due to the final bandwidth of the 
pulse alone. 



IV. CONCLUSION 

We have investigated the role of electron correlations in 
obtaining the cross section for the process of nonsequential 
two-photon double ionization of helium, using two different 
methods of approach - by projection onto uncorrected con- 
tinuum states, and by imposing a complex absorbing potential 
within the proper many-body framework. The former has the 
advantage that these states are much easier to construct than 
the true two-electron continuum states, whereas no such states 
are needed in the latter method. Moreover, within this method 
single and double ionization yields may be obtained using a 
numerical grid considerably smaller than what is necessary to 
contain the ionized part of the two-particle wave function. The 
predictions of the methods were in agreement, thus suggesting 
that projecting the ionized part of the full wave function onto 
uncorrelated continuum states does provide correct ionization 
yields in the asymptotic limit. 

Finally, we presented evidence that the cross section for this 
process remains well defined for photon energies arbitrarily 
close to the threshold at hco = 54.4 eV (from below). The 
cross section is seen to have a sharp increase in this limit. 



[1] H. Mashiko, A. Suda, and K. Midorikawa, Opt. Lett. 29, 1927 
(2004). 

[2] N. Miyamoto, M. Kamei, D. Yoshitomi, T. Kanai, T. Sekikawa, 
T. Nakajima, and S. Watanabe, Phys. Rev. Lett. 93, 083903 
(2004). 

[3] Y. Nabekawa, H. Hasegawa, E. J. Takahashi, and K. Mi- 
dorikawa, Phys. Rev. Lett. 94, 043001 (2005). 

[4] K. Hoshina, A. Hishikawa, K. Kato, T. Sako, K. Yamanouchi, 
E. J. Takahashi, Y. Nabekawa, and K. Midorikawa, J. Phys. B 
39, 813 (2006). 

[5] W. Ackermann, Nat. Photon. 1, 336 (2007). 



[6] T. Shintake et al., Nat. Photon. 2, 555 (2008). 

[7] A. A. Sorokin, S. V. Bobashev, T. Feigl, K. Tiedtke, H. Wabnitz, 

and M. Richter, Phys. Rev. Lett. 99, 213002 (2007). 
[8] R. Moshammer, Y. H. Jiang, L. Foucar, A. Rudenko, T. Ergler, 
C. D. Schroter, S. Liidemann, K. Zrost, D. Fischer, J. Titze, 
et al., Phys. Rev. Lett. 98, 203001 (2007). 
[9] T. Laarmann et al., Phys. Rev A 72, 023409 (2005). 
[10] H. Wabnitz et al., Nature 420, 482 (2002). 
[11] J. S. Parker et al., Phys. Rev. Lett. 96, 133001 (2006). 
[12] J. S. Parker, L. R. Moore, K. J. Meharg, D. Dundas, and K. T. 
Taylor, J. Phys. B 34, L69 (2001). 



8 



[13] L. R. Moore et al., J. Phys. Conf. Ser. 194, 032055 (2009). 
[14] J. Colgan and M. S. Pindzola, Phys. Rev. Lett. 88, 173002 
(2002). 

[15] L. Feng and H. W. van der Hart, J. Phys. B 36, LI (2003). 
[16] S. Laulan and H. Bachau, Phys. Rev. A 68, 013409 (2003). 
[17] B. Piraux, J. Bauer, S. Laulan, and H. Bachau, Eur. Phys. J. D 
26, 7 (2003). 

[18] S. X. Hu, J. Colgan, and L. A. Collins, J. Phys. B 38, L35 
(2005). 

[19] E. Foumouo, G. L. Kamta, G. Edah, and B. Piraux, Phys. Rev. 

A 74, 063409 (2006). 
[20] R. Shakeshaft, Phys. Rev A 76, 063405 (2007). 
[21] I. A. Ivanov and A. S. Kheifets, Phys. Rev. A 75, 033411 

(2007) . 

[22] D. A. Homer, F. Morales, T. N. Rescigno, F. Martin, and C. W. 

McCurdy, Phys. Rev. A 76, 030701(R) (2007). 
[23] L. A. A. Nikolopoulos and P. Lambropoulos, J. Phys. B 40, 

1347 (2007). 

[24] J. Feist, S. Nagele, R. Pazourek, E. Persson, B. I. Schneider, 
L. A. Collins, and J. Burgdorfer, Phys. Rev A 77, 043420 

(2008) . 

[25] X. Guan, K. Bartschat, and B. I. Schneider, Phys. Rev. A 77, 
043421 (2008). 

[26] E. Foumouo, P. Antoine, B. Piraux, L. Malegat, H. Bachau, and 

R. Shakeshaft, J. Phys. B 41, 051001 (2008). 
[27] A. Palacios, T. N. Rescigno, and C. W. McCurdy, Phys. Rev A 

79, 033402 (2009). 
[28] R. Nepstad, T. Birkeland, and M. F0rre, Phys. Rev. A 81, 

063402 (2010). 

[29] H. Hasegawa, E. J. Takahashi, Y. Nabekawa, K. L. Ishikawa, 

and K. Midorikawa, Phys. Rev. A 71, 023407 (2005). 
[30] A. A. Sorokin, M. Wellhofer, S. V. Bobashev, K. Tiedtke, and 

M. Richter, Phys. Rev. A 75, 051402(R) (2007). 
[31] A. Rudenko, L. Foucar, M. Kurka, T. Ergler, K. U. Kuhnel, 

Y. H. Jiang, A. Voitkiv, B. Najjari, A. Kheifets, S. Ludemann, 

et al., Phys. Rev. Lett. 101, 073003 (2008). 
[32] J. S. Briggs and V. Schmidt, J. Phys. B 33, Rl (2000). 
[33] L. Avaldi and A. Huetz, J. Phys. B 38, S861 (2005). 
[34] J. A. R. Samson, W. C. Stolte, Z.-X. He, J. N. Cutler, Y. Lu, and 



R. J. Bartlett, Phys. Rev A 57, 1906 (1998). 
[35] S. Selst0, and S. Kvaal, J. Phys. B 43, 065004 (2010). 
[36] K. C. Kulander, Phys. Rev. A 35, 445 (1987). 
[37] T. Fevens and H. Jiang, SIAM J. Sci. Comput. 21, 255 (1999). 
[38] B. Feuerstein and U. Thumm, J. Phys. B 36, 707 (2003). 
[39] M. Lein, E. K. U. Gross, and V. Engel, Phys. Rev Lett. 85, 4707 

(2000). 

[40] D. Bauer and F Ceccherini, Phys. Rev. A 60, 2301 (1999). 
[41] D. G. Lappas and R. van Leeuwen, J. Phys. B 31, L249 (1998). 
[42] S. L. Haan, P S. Wheeler, R. Panfili, and J. H. Eberly, Phys. 

Rev. A 66, 061402 (2002). 
[43] Z. Zhang, L.-Y. Peng, Q. Gong, and T. Morishita, Opt. Express 

18, 8976 (2010). 

[44] D. A. Horner, T. N. Rescigno, and C. W. McCurdy, Phys. Rev 

A 77, 030703(R) (2008). 
[45] P. Lambropoulos, L. A. A. Nikolopoulos, M. G. Maki'is, and 

A. Mihelic, Phys. Rev. A 78, 055402 (2008). 
[46] D. A. Horner, T. N. Rescigno, and C. W. McCurdy, Phys. Rev. 

A 81, 023410 (2010). 
[47] B. Piraux, E. Foumouo, P. Antoine, and H. Bachau, J. Phys. 

Conf. Ser. 141, 012013 (2008). 
[48] T. Birkeland, R. Nepstad, and M. F0rre, Phys. Rev. Lett. 104, 

163002 (2010). 

[49] W. S. Wilson and R. B. Lindsay, Phys. Rev. 47, 681 (1935). 
[50] V. Gorini, A. Kossakowski, and E. Sudarshan, J. Math. Phys. 

17, 821 (1976). 
[51] G. Lindblad, Comm. Math. Phys. 48, 119 (1976). 
[52] L. B. Madsen, L. A. A. Nikolopoulos, and P. Lambropoulos, 

Eur. Phys. J. D 10, 67 (2000). 
[53] E. S. Smyth, J. S. Parker, and K. Taylor, Comput. Phys. Com- 

mun. 114, 1 (1998). 
[54] H. A. Leth, L. B. Madsen, and K. M0lmer, Phys. Rev. Lett. 103, 

183601 (2009). 

[55] M. Feit, J. Fleck, and A. Steiger, J. Comput. Phys. 47, 412 
(1982). 

[56] T. Birkeland and R. Nepstad, Pyprop, 

http : / /pyprop . googlecode . com 
[57] D. E. Manolopoulos, J. Chem. Phys. 117, 9552 (2002). 



